
clear all;

% January 1, 2010 00:00:00 UTC
%tGPS = 946339215;	
%gmst = 23017.062492823; % Earth's GMSTrad at this time

% Vernal equinox in 2010: Mar 20 2010 17:32 UTC
tGPS = 953141535;
gmst = 23513.095548076; % Earth's GMSTrad at this time  

gmstRad = gmst*pi/43200;
gmstRad = 0;

%========================================================================
%--LHO detector
H1.long		= -2.08405676917;
H1.lat		= 0.81079526383;
H1.elevation	= 142.554;
H1.Ax		= -6.195e-4;
H1.PsiX		= 324.0006*pi/180.0;
H1.Ay		= 1.25e-5;
H1.PsiY		= 234.0006*pi/180.0;
H1.lmst 	= mod((gmstRad + H1.long), (2*pi));
H1.xArm		= getArm(H1.PsiX, H1.Ax, H1.long, H1.lat);
H1.yArm		= getArm(H1.PsiY, H1.Ay, H1.long, H1.lat);
H1.detTensor	= 0.5*(H1.xArm'*H1.xArm - H1.yArm'*H1.yArm);
[H1.ah, H1.bh] = CalculateDetPosDepCoeffs(H1.PsiX, H1.PsiY, H1.lat);

%-- LLO Detector
L1.long		= -1.58430937078;
L1.lat		= 0.53342313506;
L1.elevation	= -6.574;
L1.Ax		= -3.121e-4;
L1.PsiX		= 252.2835*pi/180.0;
L1.Ay		= -6.107e-4;
L1.PsiY		= 162.2835*pi/180.0;
L1.lmst 	= mod((gmstRad + L1.long), (2*pi));
L1.xArm		= getArm(L1.PsiX, L1.Ax, L1.long, L1.lat);
L1.yArm		= getArm(L1.PsiY, L1.Ay, L1.long, L1.lat);
L1.detTensor	= 0.5*(L1.xArm'*L1.xArm - L1.yArm'*L1.yArm);
[L1.ah, L1.bh] = CalculateDetPosDepCoeffs(L1.PsiX, L1.PsiY, L1.lat);

% Virgo Detector
V1.long		= 0.18333805213;
V1.lat		= 0.76151183984;
V1.elevation	= 51.884;
V1.Ax		= 0;
V1.PsiX		= 19.4326*pi/180.0;
V1.Ay		= 0;
V1.PsiY		= 289.4326*pi/180.0;
V1.lmst 	= mod((gmstRad + V1.long), (2*pi));
V1.xArm		= getArm(V1.PsiX, V1.Ax, V1.long, V1.lat);
V1.yArm		= getArm(V1.PsiY, V1.Ay, V1.long, V1.lat);
V1.detTensor	= 0.5*(V1.xArm'*V1.xArm - V1.yArm'*V1.yArm);
[V1.ah, V1.bh] = CalculateDetPosDepCoeffs(V1.PsiX, V1.PsiY, V1.lat);

% TAMA Detector
T1.long		= [139 32 9.8]*[1; 1/60; 1/3600]*pi/180; 
T1.lat		= [35 40 35.6]*[1; 1/60; 1/3600]*pi/180;
T1.elevation	= 90.0;
T1.Ax		= 0;
T1.PsiX		= 270*pi/180.0;
T1.Ay		= 0;
T1.PsiY		= 180*pi/180.0;
T1.lmst 	= mod((gmstRad + T1.long), (2*pi));
T1.xArm		= getArm(T1.PsiX, T1.Ax, T1.long, T1.lat);
T1.yArm		= getArm(T1.PsiY, T1.Ay, T1.long, T1.lat);
T1.detTensor	= 0.5*(T1.xArm'*T1.xArm - T1.yArm'*T1.yArm);
[T1.ah, T1.bh] = CalculateDetPosDepCoeffs(T1.PsiX, T1.PsiY, T1.lat);

% IndIGO Detector
I1.long		= [77 38 00]*[1; 1/60; 1/3600]*pi/180; % Bangalore
I1.lat		= [12 58 00]*[1; 1/60; 1/3600]*pi/180; % Bangalore
I1.elevation	= 920.0;
I1.Ax		= 0;
I1.PsiX		= 0;
I1.Ay		= 0;
I1.PsiY		= mod(I1.PsiX + pi/2, 2*pi);
I1.lmst 	= mod((gmstRad + I1.long), (2*pi));
I1.xArm		= getArm(I1.PsiX, I1.Ax, I1.long, I1.lat);
I1.yArm		= getArm(I1.PsiY, I1.Ay, I1.long, I1.lat);
I1.detTensor	= 0.5*(I1.xArm'*I1.xArm - I1.yArm'*I1.yArm);
[I1.ah, I1.bh] = CalculateDetPosDepCoeffs(I1.PsiX, I1.PsiY, I1.lat);
%========================================================================

% Draw up a list of n sky positions for the sources
angResolution	= .05;
skyPos		= sinusoidalMap(angResolution);
nSky		= size(skyPos,1); 

sourcePsi   = pi/4;
sourceDelta = pi/2 - skyPos(:,1);
sourceAlpha = skyPos(:,2);

[at, bt] = CalculateAtBt(H1, sourceAlpha, sourceDelta);
Fp.H1 = (at * cos(2*sourcePsi)) + (bt * sin(2*sourcePsi));
Fc.H1 = (bt * cos(2*sourcePsi)) - (at * sin(2*sourcePsi));

clear at bt
[at, bt] = CalculateAtBt(V1, sourceAlpha, sourceDelta);
Fp.V1 = (at * cos(2*sourcePsi)) + (bt * sin(2*sourcePsi));
Fc.V1 = (bt * cos(2*sourcePsi)) - (at * sin(2*sourcePsi));

clear at bt
[at, bt] = CalculateAtBt(T1, sourceAlpha, sourceDelta);
Fp.T1 = (at * cos(2*sourcePsi)) + (bt * sin(2*sourcePsi));
Fc.T1 = (bt * cos(2*sourcePsi)) - (at * sin(2*sourcePsi));

clear at bt
[at, bt] = CalculateAtBt(I1, sourceAlpha, sourceDelta);
Fp.I1 = (at * cos(2*sourcePsi)) + (bt * sin(2*sourcePsi));
Fc.I1 = (bt * cos(2*sourcePsi)) - (at * sin(2*sourcePsi));

%-- Plotting part
projectionType = 'miller';
xproject( skyPos, sqrt(Fp.I1.^2+Fc.I1.^2), projectionType ); 

% load the worldmap coastlines:
mapstruct = defaultm(projectionType);
mapstruct.angleunits = 'degrees';
mapstruct.origin = [0 0 0];
mapstruct.falseeasting = 0;
mapstruct.falsenorthing = 0;
mapstruct.scalefactor = 1;
load coast;
[xc, yc] = mfwdtran(mapstruct, lat, long);
hold on;
plot(xc, yc, 'k-', 'LineWidth', 1); % coastlines
hold off

% Plot the location of all the detectors
detLat = [H1.lat L1.lat V1.lat I1.lat T1.lat]*180/pi;
detLon = [H1.long L1.long V1.long I1.long T1.long]*180/pi;
[detx, dety] = mfwdtran(mapstruct, detLat, detLon);
hold on;
plot(detx, dety, 'mo', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
hold off;

box on;
set(gca, 'FontSize', 12);
xlabel('Longitude / radian', 'FontSize', 14);
ylabel('Latitude / radian', 'FontSize', 14);
caxis([0 1.6])
























